home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / zlatrd.f < prev    next >
Text File  |  1997-06-25  |  11KB  |  281 lines
  1.       SUBROUTINE ZLATRD( UPLO, N, NB, A, LDA, E, TAU, W, LDW )
  2. *
  3. *  -- LAPACK auxiliary routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     September 30, 1994
  7. *
  8. *     .. Scalar Arguments ..
  9.       CHARACTER          UPLO
  10.       INTEGER            LDA, LDW, N, NB
  11. *     ..
  12. *     .. Array Arguments ..
  13.       DOUBLE PRECISION   E( * )
  14.       COMPLEX*16         A( LDA, * ), TAU( * ), W( LDW, * )
  15. *     ..
  16. *
  17. *  Purpose
  18. *  =======
  19. *
  20. *  ZLATRD reduces NB rows and columns of a complex Hermitian matrix A to
  21. *  Hermitian tridiagonal form by a unitary similarity
  22. *  transformation Q' * A * Q, and returns the matrices V and W which are
  23. *  needed to apply the transformation to the unreduced part of A.
  24. *
  25. *  If UPLO = 'U', ZLATRD reduces the last NB rows and columns of a
  26. *  matrix, of which the upper triangle is supplied;
  27. *  if UPLO = 'L', ZLATRD reduces the first NB rows and columns of a
  28. *  matrix, of which the lower triangle is supplied.
  29. *
  30. *  This is an auxiliary routine called by ZHETRD.
  31. *
  32. *  Arguments
  33. *  =========
  34. *
  35. *  UPLO    (input) CHARACTER
  36. *          Specifies whether the upper or lower triangular part of the
  37. *          Hermitian matrix A is stored:
  38. *          = 'U': Upper triangular
  39. *          = 'L': Lower triangular
  40. *
  41. *  N       (input) INTEGER
  42. *          The order of the matrix A.
  43. *
  44. *  NB      (input) INTEGER
  45. *          The number of rows and columns to be reduced.
  46. *
  47. *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
  48. *          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
  49. *          n-by-n upper triangular part of A contains the upper
  50. *          triangular part of the matrix A, and the strictly lower
  51. *          triangular part of A is not referenced.  If UPLO = 'L', the
  52. *          leading n-by-n lower triangular part of A contains the lower
  53. *          triangular part of the matrix A, and the strictly upper
  54. *          triangular part of A is not referenced.
  55. *          On exit:
  56. *          if UPLO = 'U', the last NB columns have been reduced to
  57. *            tridiagonal form, with the diagonal elements overwriting
  58. *            the diagonal elements of A; the elements above the diagonal
  59. *            with the array TAU, represent the unitary matrix Q as a
  60. *            product of elementary reflectors;
  61. *          if UPLO = 'L', the first NB columns have been reduced to
  62. *            tridiagonal form, with the diagonal elements overwriting
  63. *            the diagonal elements of A; the elements below the diagonal
  64. *            with the array TAU, represent the  unitary matrix Q as a
  65. *            product of elementary reflectors.
  66. *          See Further Details.
  67. *
  68. *  LDA     (input) INTEGER
  69. *          The leading dimension of the array A.  LDA >= max(1,N).
  70. *
  71. *  E       (output) DOUBLE PRECISION array, dimension (N-1)
  72. *          If UPLO = 'U', E(n-nb:n-1) contains the superdiagonal
  73. *          elements of the last NB columns of the reduced matrix;
  74. *          if UPLO = 'L', E(1:nb) contains the subdiagonal elements of
  75. *          the first NB columns of the reduced matrix.
  76. *
  77. *  TAU     (output) COMPLEX*16 array, dimension (N-1)
  78. *          The scalar factors of the elementary reflectors, stored in
  79. *          TAU(n-nb:n-1) if UPLO = 'U', and in TAU(1:nb) if UPLO = 'L'.
  80. *          See Further Details.
  81. *
  82. *  W       (output) COMPLEX*16 array, dimension (LDW,NB)
  83. *          The n-by-nb matrix W required to update the unreduced part
  84. *          of A.
  85. *
  86. *  LDW     (input) INTEGER
  87. *          The leading dimension of the array W. LDW >= max(1,N).
  88. *
  89. *  Further Details
  90. *  ===============
  91. *
  92. *  If UPLO = 'U', the matrix Q is represented as a product of elementary
  93. *  reflectors
  94. *
  95. *     Q = H(n) H(n-1) . . . H(n-nb+1).
  96. *
  97. *  Each H(i) has the form
  98. *
  99. *     H(i) = I - tau * v * v'
  100. *
  101. *  where tau is a complex scalar, and v is a complex vector with
  102. *  v(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in A(1:i-1,i),
  103. *  and tau in TAU(i-1).
  104. *
  105. *  If UPLO = 'L', the matrix Q is represented as a product of elementary
  106. *  reflectors
  107. *
  108. *     Q = H(1) H(2) . . . H(nb).
  109. *
  110. *  Each H(i) has the form
  111. *
  112. *     H(i) = I - tau * v * v'
  113. *
  114. *  where tau is a complex scalar, and v is a complex vector with
  115. *  v(1:i) = 0 and v(i+1) = 1; v(i+1:n) is stored on exit in A(i+1:n,i),
  116. *  and tau in TAU(i).
  117. *
  118. *  The elements of the vectors v together form the n-by-nb matrix V
  119. *  which is needed, with W, to apply the transformation to the unreduced
  120. *  part of the matrix, using a Hermitian rank-2k update of the form:
  121. *  A := A - V*W' - W*V'.
  122. *
  123. *  The contents of A on exit are illustrated by the following examples
  124. *  with n = 5 and nb = 2:
  125. *
  126. *  if UPLO = 'U':                       if UPLO = 'L':
  127. *
  128. *    (  a   a   a   v4  v5 )              (  d                  )
  129. *    (      a   a   v4  v5 )              (  1   d              )
  130. *    (          a   1   v5 )              (  v1  1   a          )
  131. *    (              d   1  )              (  v1  v2  a   a      )
  132. *    (                  d  )              (  v1  v2  a   a   a  )
  133. *
  134. *  where d denotes a diagonal element of the reduced matrix, a denotes
  135. *  an element of the original matrix that is unchanged, and vi denotes
  136. *  an element of the vector defining H(i).
  137. *
  138. *  =====================================================================
  139. *
  140. *     .. Parameters ..
  141.       COMPLEX*16         ZERO, ONE, HALF
  142.       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ),
  143.      $                   ONE = ( 1.0D+0, 0.0D+0 ),
  144.      $                   HALF = ( 0.5D+0, 0.0D+0 ) )
  145. *     ..
  146. *     .. Local Scalars ..
  147.       INTEGER            I, IW
  148.       COMPLEX*16         ALPHA
  149. *     ..
  150. *     .. External Subroutines ..
  151.       EXTERNAL           ZAXPY, ZGEMV, ZHEMV, ZLACGV, ZLARFG, ZSCAL
  152. *     ..
  153. *     .. External Functions ..
  154.       LOGICAL            LSAME
  155.       COMPLEX*16         ZDOTC
  156.       EXTERNAL           LSAME, ZDOTC
  157. *     ..
  158. *     .. Intrinsic Functions ..
  159.       INTRINSIC          DBLE, MIN
  160. *     ..
  161. *     .. Executable Statements ..
  162. *
  163. *     Quick return if possible
  164. *
  165.       IF( N.LE.0 )
  166.      $   RETURN
  167. *
  168.       IF( LSAME( UPLO, 'U' ) ) THEN
  169. *
  170. *        Reduce last NB columns of upper triangle
  171. *
  172.          DO 10 I = N, N - NB + 1, -1
  173.             IW = I - N + NB
  174.             IF( I.LT.N ) THEN
  175. *
  176. *              Update A(1:i,i)
  177. *
  178.                A( I, I ) = DBLE( A( I, I ) )
  179.                CALL ZLACGV( N-I, W( I, IW+1 ), LDW )
  180.                CALL ZGEMV( 'No transpose', I, N-I, -ONE, A( 1, I+1 ),
  181.      $                     LDA, W( I, IW+1 ), LDW, ONE, A( 1, I ), 1 )
  182.                CALL ZLACGV( N-I, W( I, IW+1 ), LDW )
  183.                CALL ZLACGV( N-I, A( I, I+1 ), LDA )
  184.                CALL ZGEMV( 'No transpose', I, N-I, -ONE, W( 1, IW+1 ),
  185.      $                     LDW, A( I, I+1 ), LDA, ONE, A( 1, I ), 1 )
  186.                CALL ZLACGV( N-I, A( I, I+1 ), LDA )
  187.                A( I, I ) = DBLE( A( I, I ) )
  188.             END IF
  189.             IF( I.GT.1 ) THEN
  190. *
  191. *              Generate elementary reflector H(i) to annihilate
  192. *              A(1:i-2,i)
  193. *
  194.                ALPHA = A( I-1, I )
  195.                CALL ZLARFG( I-1, ALPHA, A( 1, I ), 1, TAU( I-1 ) )
  196.                E( I-1 ) = ALPHA
  197.                A( I-1, I ) = ONE
  198. *
  199. *              Compute W(1:i-1,i)
  200. *
  201.                CALL ZHEMV( 'Upper', I-1, ONE, A, LDA, A( 1, I ), 1,
  202.      $                     ZERO, W( 1, IW ), 1 )
  203.                IF( I.LT.N ) THEN
  204.                   CALL ZGEMV( 'Conjugate transpose', I-1, N-I, ONE,
  205.      $                        W( 1, IW+1 ), LDW, A( 1, I ), 1, ZERO,
  206.      $                        W( I+1, IW ), 1 )
  207.                   CALL ZGEMV( 'No transpose', I-1, N-I, -ONE,
  208.      $                        A( 1, I+1 ), LDA, W( I+1, IW ), 1, ONE,
  209.      $                        W( 1, IW ), 1 )
  210.                   CALL ZGEMV( 'Conjugate transpose', I-1, N-I, ONE,
  211.      $                        A( 1, I+1 ), LDA, A( 1, I ), 1, ZERO,
  212.      $                        W( I+1, IW ), 1 )
  213.                   CALL ZGEMV( 'No transpose', I-1, N-I, -ONE,
  214.      $                        W( 1, IW+1 ), LDW, W( I+1, IW ), 1, ONE,
  215.      $                        W( 1, IW ), 1 )
  216.                END IF
  217.                CALL ZSCAL( I-1, TAU( I-1 ), W( 1, IW ), 1 )
  218.                ALPHA = -HALF*TAU( I-1 )*ZDOTC( I-1, W( 1, IW ), 1,
  219.      $                 A( 1, I ), 1 )
  220.                CALL ZAXPY( I-1, ALPHA, A( 1, I ), 1, W( 1, IW ), 1 )
  221.             END IF
  222. *
  223.    10    CONTINUE
  224.       ELSE
  225. *
  226. *        Reduce first NB columns of lower triangle
  227. *
  228.          DO 20 I = 1, NB
  229. *
  230. *           Update A(i:n,i)
  231. *
  232.             A( I, I ) = DBLE( A( I, I ) )
  233.             CALL ZLACGV( I-1, W( I, 1 ), LDW )
  234.             CALL ZGEMV( 'No transpose', N-I+1, I-1, -ONE, A( I, 1 ),
  235.      $                  LDA, W( I, 1 ), LDW, ONE, A( I, I ), 1 )
  236.             CALL ZLACGV( I-1, W( I, 1 ), LDW )
  237.             CALL ZLACGV( I-1, A( I, 1 ), LDA )
  238.             CALL ZGEMV( 'No transpose', N-I+1, I-1, -ONE, W( I, 1 ),
  239.      $                  LDW, A( I, 1 ), LDA, ONE, A( I, I ), 1 )
  240.             CALL ZLACGV( I-1, A( I, 1 ), LDA )
  241.             A( I, I ) = DBLE( A( I, I ) )
  242.             IF( I.LT.N ) THEN
  243. *
  244. *              Generate elementary reflector H(i) to annihilate
  245. *              A(i+2:n,i)
  246. *
  247.                ALPHA = A( I+1, I )
  248.                CALL ZLARFG( N-I, ALPHA, A( MIN( I+2, N ), I ), 1,
  249.      $                      TAU( I ) )
  250.                E( I ) = ALPHA
  251.                A( I+1, I ) = ONE
  252. *
  253. *              Compute W(i+1:n,i)
  254. *
  255.                CALL ZHEMV( 'Lower', N-I, ONE, A( I+1, I+1 ), LDA,
  256.      $                     A( I+1, I ), 1, ZERO, W( I+1, I ), 1 )
  257.                CALL ZGEMV( 'Conjugate transpose', N-I, I-1, ONE,
  258.      $                     W( I+1, 1 ), LDW, A( I+1, I ), 1, ZERO,
  259.      $                     W( 1, I ), 1 )
  260.                CALL ZGEMV( 'No transpose', N-I, I-1, -ONE, A( I+1, 1 ),
  261.      $                     LDA, W( 1, I ), 1, ONE, W( I+1, I ), 1 )
  262.                CALL ZGEMV( 'Conjugate transpose', N-I, I-1, ONE,
  263.      $                     A( I+1, 1 ), LDA, A( I+1, I ), 1, ZERO,
  264.      $                     W( 1, I ), 1 )
  265.                CALL ZGEMV( 'No transpose', N-I, I-1, -ONE, W( I+1, 1 ),
  266.      $                     LDW, W( 1, I ), 1, ONE, W( I+1, I ), 1 )
  267.                CALL ZSCAL( N-I, TAU( I ), W( I+1, I ), 1 )
  268.                ALPHA = -HALF*TAU( I )*ZDOTC( N-I, W( I+1, I ), 1,
  269.      $                 A( I+1, I ), 1 )
  270.                CALL ZAXPY( N-I, ALPHA, A( I+1, I ), 1, W( I+1, I ), 1 )
  271.             END IF
  272. *
  273.    20    CONTINUE
  274.       END IF
  275. *
  276.       RETURN
  277. *
  278. *     End of ZLATRD
  279. *
  280.       END
  281.